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We present the baseline multimessenger analysis method for the joint observations of gravitational 
waves (GW) and high-energy neutrinos (HEN), together with a detailed analysis of the expected 
science reach of the joint search. The analysis method combines data from GW and HEN detectors, 
and uses the blue-luminosity-weighted distribution of galaxies. We derive expected GW+HEN 
source rate upper limits for a wide range of source parameters covering several emission models. 
Using published sensitivities of externally triggered searches, we derive joint upper limit estimates 
both for the ongoing analysis with the initial LIGO-Virgo GW detectors with the partial IceCube 
detector (22 strings) HEN detector and for projected results to advanced LIGO-Virgo detectors with 
the completed IceCube (86 strings). We discuss the constraints these upper limits impose on some 
existing GW+HEN emission models. 



I. INTRODUCTION 

The observation of gravitational- waves (GWs) and 
neutrinos is entering a new and promising era with newly 
constructed detectors. GW observatories such as LIGO 
[1], Virgo [2] and GEO [3] will be upgraded to second 
generation detectors within the next few years. Another 
advanced GW detector, LCGT [4], is being constructed 
in Japan, while LIGO is considering plans for a third 
observatory in India [5]. 

Both the emission mechanism and detection method of 
neutrinos can depend greatly on their energy, therefore 
it is worthwhile to consider different neutrino subgroups. 
So far, only astrophysical MeV neutrinos (in addition to 
those from the Sun) have been detected, and only in one 
case, from supernova SN 1987A [6, 7]. Here we concen- 
trate specifically on high-energy neutrinos (HENs): neu- 
trinos > 100 GeV that carry information about the par- 
ticle acceleration region of astrophysical sources [8]. 
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HEN observatories currently in operation are Ice- 
Cube [9], a cubic-kilometer detector at the geographic 
South Pole, and Antares [10] in the Mediterranean sea. 
Antares is being upgraded to a cubic-kilometer detector 
called KM3NeT in the following years [11]. A third HEN 
detector operating at the lake Baikal is also planned to 
be upgraded to a km 3 volume [12]. 

The joint (multimessenger) analysis of GW and HEN 
observations presents multiple advantages over single 
messenger analyses. Since both GWs and HENs are 
weakly interacting, their detection requires exception- 
ally high sensitivity. Search sensitivity can be greatly 
increased by combining data from GW and HEN detec- 
tors. The multimessenger approach could also signifi- 
cantly add to our understanding of the underlying mecha- 
nisms that create the astrophysical sources emitting both 
signals [13]. Aso et al [14] designed a GW+HEN multi- 
messenger search algorithm that uses the two LIGO de- 
tectors and IceCube to look for spatially and temporally 
coincident events. After requiring temporal coincidence 
between the GW and HEN signals within a given time 
window, the method reconstructs a ring on the sky based 
on the measured time delay between the signal arrival 
times in the LIGO detectors, and requires the direction 
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of the neutrino signal to overlap with this ring. Aso et at. 
show that such coincidence is extremely unlikely to arise 
from the background, making the requirement of coinci- 
dence very effective in reducing the false alarm rate of 
the joint measurement. 

Pradier [15] considered a joint search with initial Virgo 
and Antares, and discussed the feasibility of a time 
coincident search using these two detectors. Pradier 
discussed microquasars and flares from soft gamma re- 
peaters as plausible galactic GW+HEN sources. He also 
considered the detectability of quantum gravity effects 
by measuring the time delay between the arrival of GW 
and HEN signals. 

We are close to the milestone of finishing the first co- 
incident search for GWs and HENs for the initial LIGO- 
Virgo (S5/VSR1 science runs) and the partial Antares 
detector in its 5-string configuration. The analysis, a 
simpler version of what is presented below, uses the di- 
rectional distribution and the time of arrival of HENs to 
trigger a GW follow-up analysis, similar to the analy- 
sis used for GW follow-up searches of GRBs (e.g. [16]). 
There are ~200 neutrino triggers from the Antares, 
most of which are detected by digital optical modules 
(DOMs) two strings, while 13 neutrinos are detected by 
DOMs on three strings. The first scientific results of this 
search will be published soon. 

In this article we introduce a joint GW and HEN analy- 
sis algorithm that combines the observations of a network 
of GW detectors and a HEN detector. Besides looking 
for astrophysical GW+HEN messengers, the search al- 
gorithm can also be used to derive upper limits on the 
population and flux of GW+HEN sources. We estimate 
the anticipated science reach for initial and advanced de- 
tectors, and discuss some of the existing emission models 
and how they can be constrained in the event of non- 
detection. 

The distribution of astrophysical GW+HEN sources at 
detectable distances is not uniform. This can be used in a 
joint search algorithm to reduce false coincidences and in- 
crease sensitivity. One can weigh event candidates based 
on the expected source density in their direction. The 
method utilizes the blue-luminosity-weighted galaxy dis- 
tribution to favor astrophysical sources over background 
coincidences. Blue luminosity is a good tracer of the rate 
of star formation and therefore CCSN rate (e.g. [17]). 
Long GRB rates in galaxies also correlate with the galax- 
ies' blue luminosity [18] (typically long GRBs are more 
likely to occur in smaller, bluer galaxies compared to CC- 
SNe. The method weighs source directions with the ex- 
pected detectable source rate from the given directions, 
assuming that source distribution follows blue-luminosity 
distribution. We take the blue-luminosity distribution 
of galaxies up to 40 Mpc from the Gravitational Wave 
Galaxy Catalogue (GWGC) [19]. 

The article is organized as follows. In section II we dis- 
cuss some anticipated astrophysical sources of GWs and 
HENs. Section III discusses the expected science reach 
of joint GW+HEN searches by presenting the interpreta- 



tion of expected results for constraining existing emission 
and population models. Section IV describes GW and 
HEN detectors and data. In section V we introduce the 
baseline joint GW+HEN analysis, also discussing its re- 
lation to the presented science reach. Section VI presents 
a summary of the method and the science reach. 



II. ASTROPHYSICAL SOURCES 

GWs and HENs may originate from a number of com- 
mon sources. Plausible sources include gamma-ray bursts 
(GRBs) [20-28], core-collapse supernovae (CCSNe), soft 
gamma repeaters [29-31] and microquasars [13, 32]. For 
a joint GW+HEN search, potentially the most interest- 
ing sources are those which are difficult to detect using 
electromagnetic (EM) telescopes. Below we concentrate 
on GRBs as plausible GW+HEN sources. We discuss the 
scenarios in which GRBs have limited EM emission. 

GRBs are exceptionally luminous gamma-ray flashes 
of cosmic origin [13]. They are thought to originate from 
at least two distinct types of astrophysical sources. The 
core-collapse of massive stars is thought to be the pro- 
genitor of at least some long-soft GRBs [33] , while short- 
hard GRBs are usually associated with the mergers of 
compact binaries, such as two neutron stars (NS), or NS 
- black hole (BH) binaries [34]. 

Both types of GRB progenitors are expected to emit 
GWs. In the collapsar scenario, transient GW bursts are 
likely to be emitted by a variety of processes inside the 
star. These processes include rotating core collapse and 
bounce, non-axisymmetric rotational instabilities, post- 
bounce convective overturn, or non-radial PNS pulsa- 
tions [35]. Other mechanisms, such as the fragmenta- 
tion of collapsar accretion disks [36-38], or suspended 
accretion [39], might also play an important role in GW 
emission. 

Compact binary mergers are expected to be strong 
GW emitters in the sensitive frequency band of Earth- 
based GW detectors [40]. Binary systems lose angular 
momentum due to the emission of GWs. In the inspiral 
phase, the distance between the two compact objects de- 
creases, while the rotational frequency increases. The 
process is expected to continue until the two objects 
merge, an event which is anticipated to involve an in- 
termediate stage of a central quasi- axi-symmetric object 
surrounded by an accretion disk [34] , producing a strong 
GW transient [40]. 

Predictions on the energy radiated away through GWs 
from GRB progenitors vary over several orders of magni- 
tude, depending on the emission model considered. Nu- 
merical simulations of various non-rotating CCSN pro- 
genitor models give GW emission of 10 -8 — 10 -4 M c 2 
[35, 41], at typical frequencies of 500 — 700 Hz. Re- 
cent simulations of rotating CCSN progenitors with ini- 
tial conditions chosen to resemble likely long- GRB pro- 
genitors show GW emission of E GW ~ 10 _T M c 2 [42], 
with characteristic frequencies of f c ~ 500 — 1000 Hz. 
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Such emissions are detectable with advanced GW detec- 
tors from sources within our Galaxy. Other, analytical 
and numerical models of long-GRBs that consider GW 
production by matter accretion around a central black 
hole tend to predict significantly stronger emission of 
GWs. For instance the fragmentation of collapsar accre- 
tion disks [36-38] could emit 1(T 3 - 1(T 2 M c 2 in GWs, 
potentially in the most sensitive frequency band of LIGO- 
Virgo (~ 150 Hz). Recent simulations of black hole-torus 
systems show that non-axisymmetric instabilities in such 
systems may produce strong, quasi-periodic GW signals 
in the sensitive frequency band of ground based detectors, 
potentially detectable with advanced detectors from up 
to ~ 100 Mpc [43, 44]. According to the suspended ac- 
cretion model [39], GW emission up to 10 -2 — 10 _1 M c 2 
is possible in the most sensitive frequency band. Mod- 
els for the likely progenitors of short GRBs, i.e. black 
hole - neutron star (NS) or NS-NS binaries, are expected 
to emit up to 10 _1 M c 2 in GWs in the most sensitive 
frequency band [45]. 

HENs (> 100 GeV neutrinos) are thought to be pro- 
duced by various sources, including GRBs, as described 
by the internal shock model (e.g. [20-22, 25]). According 
to the model, a central engine accelerates protons and 
electrons to relativistic velocities through Fermi acceler- 
ation [20, 46, 47]. Relativistic electrons emit gamma-rays 
through synchrotron radiation, while relativistic protons 
interact with these gamma-rays (jp) or with other pro- 
tons, (pp) producing charged pions (7r ± ) and kaons (i^ ± ). 
Charged pions and kaons create HENs through the decay 
process [48] 

7r ± ,Jf ± -^^ + 1/^) (1) 

Muons may further decay and emit an additional muon- 
neutrino and electron-neutrino. However, muons may 
interact or undergo radiative cooling before they decay, 
in which case the contribution of secondary high-energy 
neutrinos will be reduced [48, 49]. 

To characterize the HEN flux of astrophysical sources, 
we use the average number of detected HENs from a HEN 
source at 10 Mpc in a random direction, denoted by n REN . 
This number depends on the sensitivity of the neutrino 
detector. 

The Waxman-Bahcall model [20], the benchmark 
model of HEN emission from GRBs, predicts about 
n HE N ~ 100 neutrinos detected in a km 3 detector for 
a typical GRB at 10 Mpc (e.g. [50]). Another model 
of HEN emission from mildly relativistic jets of core- 
collapse supernovae, and potentially from choked GRBs 
(below), predicts HEN emission of n HE N ~ 10 [48] (note 
that the result presented in [48] is three times higher, as 
it does not take into account neutrino flavor mixing) . Ho- 
riuchi and Ando [51] estimate n UEN from reverse shocks in 
mildly relativistic jets to be n HE N ~ 0.7 — 7 for a km 3 neu- 
trino detector (after taking into account neutrino flavor 
mixing). Meszaros and Waxman [22] predict a emission 
of 10 _1 - 10 detected neutrinos by a km 3 neutrino de- 
tector at a cosmological distance of z ~ 1 for collapsars. 



Razzaque et al [49] obtain n HE N ~ 0.15 for supernovae 
with mildly relativistic jets with jet energy of E ~ 10 51,5 
erg. Razzaque et al. find this result to scale linearly with 
jet energy (i.e. it is ~ xlO higher for typical hypernovae). 

The joint search for astrophysical GW+HEN sources 
is of special importance for sources which are barely or 
not at all observable through other messengers. Below 
we discuss some scenarios in which EM emission from 
GRBs is limited. The discovery of such sources is one 
of the most valuable scientific goals of joint GW+HEN 
searches. 



A. Choked GRBs 

Core-collapse supernovae (CCSN) can have observable 
gamma-ray emission only if the relativistic outflow from 
the central engine, that is responsible for the production 
of gamma-rays, breaks out of the star [51]. The outflow 
can advance only as long as it is driven by the central 
engine. If the duration of the activity of the central en- 
gine is shorter than the breakout time of the outflow, the 
outflow is choked, resulting in a choked GRB [52]. 

HENs, due to their weak interaction with matter, 
can escape from inside the stellar envelope, from depths 
gamma-rays cannot. Consequently choked GRBs, simi- 
larly to "successful" GRBs, may be significant sources of 
HENs [51, 53]. Choked GRBs are expected to emit GWs 
similarly to successful GRBs. 



B. Low-luminosity GRBs 

Low- luminosity (LL) GRBs [54-61] form a sub-class of 
long GRBs with sub-energetic gamma-ray emission. Of 
the six detected long-GRBs with spectroscopically con- 
firmed SN associations, four are of the LL-GRB variety 
[62]. Although in general both of these classes of GRBs 
have been associated with luminous Type Ic SNe, the 
less energetic SN that accompanied low-luminous burst 
GRB 060218 suggests that the connection may extend 
towards lower energy SN explosions (e.g. [58]). In addi- 
tion, the relative close proximity of observed LL-GRBs 
implies a much higher rate of occurrence than that of 
canonical high-luminosity (HL) GRBs [58, 60]. Due to 
this higher population rate, the total diffuse HEN flux 
from LL-GRBs can be comparable or even surpass the 
flux from conventional HL GRBs [28, 63, 64]. While 
individual LL GRBs are less luminous in neutrinos, the 
higher population rate makes LL-GRBs valuable sources 
for GW+HEN searches. 



C. Unknown Sources and Mechanisms 

A potential advantage of the search for astrophysical 
GW and HEN signals is the discovery of previously unan- 
ticipated sources or mechanisms. Such mechanisms in- 
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elude, for example, HEN emission at a larger beaming 
angle than the observed beaming angle for gamma rays. 

D. Source population 

Some of the most interesting GW+HEN sources, CC- 
SNe (SNe Ib/c) are expected to be similarly or somewhat 
less abundant [48] than galactic SNe, with an estimated 
rate of ~1 /century /Milky Way equivalent (MWE) galaxy 
[65]. Additionally, late-time radio observations of super- 
novae indicate that < 1% of SNe have mildly relativistic 
jets [66]. The expected rate of another interesting pop- 
ulation, LL-GRBs, is - 3 x 1(T 5 yr" 1 / galaxy (-300 
Gpc _3 yr _1 ) [67], although this value is highly uncertain 
[68]. The rate of HL-GRBs is several orders of magni- 
tudes smaller (~ 1 Gpc _3 yr _1 ) [67]. 

E. Time Delay between Gravitational- waves and 
High-energy Neutrinos 

GWs and HENs are though to be emitted through dif- 
ferent emission mechanisms by a joint source. Due to 
this difference the two messengers will likely be observed 
with a time delay. Further time delay can arise from the 
randomness of HEN-detection throughout the emission 
period. The time difference between the observation of 
GWs and HENs can be an important factor in the inter- 
pretation of multimessenger signals [69] . For instance an 
HEN detected prior to a GW signal from a GRB would 
indicate precursor activity in the GRB, while the time 
delay of the earliest HEN after the GW signature of a 
collapsar may be indicative of jet propagation within the 
stellar envelope. Besides interpretation, an upper bound 
on the temporal difference between the observation of 
GWs and HENs is an important parameter in designing 
a joint search algorithm (see Section V A below). 

Baret et al. [69] recently published an estimate on the 
upper bound of the time delay between GWs and HENs 
from GRBs. Their analysis focused on GRBs as arguably 
the most promising multimessenger GW+HEN sources. 
They obtained a conservative ~ 500 s time window that 
took into account processes motivated by current GRB 
models. The duration of each process was determined 
based on electromagnetic observational data. 



III. SCIENCE REACH 

In this section we investigate the constraints one can 
introduce with the GW+HEN search on the population 
of astrophysical GW+HEN sources. Below we first es- 
timate the expected population upper limits from the 
GW+HEN search as a function of source parameters, 
after which we interpret these constraints. The science 
reach analysis presented here follows the method of Bar- 
tos et al. [93], that we outline below. 



In determining the GW+HEN population upper limit 
we assume standard GW+HEN sources with the same in- 
trinsic emission. Limits based on a fixed average bright- 
ness are conservative compared to those using any other 
brightness distribution. We consider maximum one HEN 
detected for each source. This is a reasonable (and con- 
servative) assumption given that there has been no dis- 
covery with neutrino detectors. We introduce the exclu- 
sion distance D^™% EN , which is the maximum distance 
that satisfies the following criterion: for an astrophysi- 
cal GW+HEN burst at a distance < D^ EN in a typical 
direction and with one detected HEN, the probability 
that it is more significant (Eq. 17) than the loudest ob- 
served event of the GW+HEN search is > 50%. This dis- 
tance depends on the total (isotropic-equivalent) energy 
emitted in GWs (E^) of a GW+HEN source. Using 
this distance, we calculate the minimum astrophysical 
GW+HEN source rate (i.e. population) that would have 
produced at least one detected astrophysical HEN signal 
with > 90% probability. This source rate will be the rate 
upper limit. 

To estimate the expected results with the GW+HEN 
search, we approximate the sensitivity of the GW+HEN 
algorithm with that of published externally triggered GW 
searches (e.g. [16, 70]). This is a reasonable (and con- 
servative) approximation if one chooses the threshold for 
GW and HEN trigger selection such that there will be 
O(l) spatially and temporally coincident GW and HEN 
signals for the duration of the measurement. In an ex- 
ternally triggered search for GW bursts in coincidence 
with GRBs, Abbott et al. [16] obtained a median upper 
bound of h^sl ~ 3.8 x 10 -22 strain root sum square using 
the initial LIGO- Virgo GW detector network (S5/VSR1 
science run). They used a sine-Gaussian waveform with 
characteristic frequency of ~ 150 Hz, which is in the 
most sensitive frequency band of the GW detectors. This 
upper limit corresponds to the minimal GW signal am- 
plitude that, with 90% confidence, produces larger sig- 
nificance than the loudest joint GW+HEN event in the 
real data measured in coincidence with an external trig- 
ger. To estimate the performance of advanced detectors 
(advanced LIGO-Virgo), we estimate their median strain 
upper bound as 0.1 times that of initial detectors (i.e. 
3.8 x 10 -23 ). We note here that with additional advanced 
detectors, such as LIGO- Australia [71] and LCGT [4], 
the sensitivity of the GW detector network will further 
increase. For comparison, we note that the upper bound 
obtained with the all-sky GW search of Abadie et al. 
[72] for sine-Gaussian signals at ~ 150 Hz with the ini- 
tial LIGO-Virgo detector network is h^~ sky « 6 x 10~ 22 . 
The all-sky search corresponds to the lower limit for the 
sensitivity of GW searches as no additional information 
is used besides GW data. 

For the GW+HEN search we introduce the upper 
bound ^^ HEN , and estimate this upper bound to be 
h?7s H ™ ~ Kfs- We assume that a GW signal with 
h rss > fo^ HEN in coincidence with a detected astrophysi- 
cal HEN would be detected by the joint GW+HEN search 
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FIG. 1. Expected GW+HEN source population upper limits for IceCube-22 coincident with initial LIGO- Virgo (left) and 
IceCube-86 coincident with advanced the LIGO- Virgo detectors (right; courtesy of [93]), with one year of coincident measure- 
ment time. The results take into account the blue-luminosity-weighted galaxy distribution. The x-axis represents the GW 
energy output of a standard source. The y-axis represents the number of detected neutrinos from a standard source at 10 Mpc. 
The color scale shows the obtained source rate upper limit R UL in logarithmic units of number of sources per (Milky Way 
equivalent) galaxy per year. On both plots, the two horizontal lines (scaled for detector sensitivity) show the Waxman-Bahcall 
emission model [50] (higher) and the HEN emission model of Ando and Beacom [48] for reverse shocks in mildly relativistic 
supernova jets / choked GRBs (lower). 



with > 90% probability. Given h™ HEN (i.e. the ampli- 
tude at the detector) and Eq^ (i.e. the amplitude at the 
source), we can calculate the radius within which there 
was likely (P > 90%) no GW+HEN event from which 
a HEN was detected. This distance, averaged over all 
directions on the sky, is [16] 



D G 



12 Mpc 



zpiso 



10- 2 M C 2 



1/2 



hi 



h GWHEN 
IL rss 



(2) 



From the fact that no astrophysical HEN was de- 
tected from a GW+HEN source within D GWHEN , we ob- 
tain a source rate upper limit as the highest source rate 
that would have produced at least one detected neutrino 
within D GWHEN with < 90% probability. Assuming a Pois- 
sonian source number distribution, this corresponds to 
an average source rate of 2.3 over one year long measure- 
ment. We denote this source rate upper limit by i? UL . To 
obtain i? UL , we calculate the average number of sources 
within £) GWHEN over the duration of the measurement 
from which at least one neutrino has been detected. i? UL 
will be the rate that corresponds to 2.3 detected sources 
on average. The rate depends on the blue-luminosity- 
weighted galaxy distribution within £) GWHEN (see section 
IV D), as well as the neutrino flux (n UEN ) from a standard 
source. 

To scale theoretical expectations on the HEN rate for 
km 3 detectors to expectations for the IceCube 22 string 
(hereafter IceCube-22) detector (which was operating the 
same time as the initial LIGO- Virgo detectors (S5/VSR1 
science run)), we estimate IceCube-22 to be approxi- 



mately 10 times less sensitive. 

The estimated source rate upper limit is dependent on 
the beaming of HEN emission (beaming is less significant 
for GWs). The beaming of HEN emission is uncertain, 
but it is probably similar to the beaming of gamma rays, 
as the two emission mechanisms are connected. For this 
reason we use the gamma-ray beaming factor obtained 
for LL-GRBs, estimated to be less than 14 [67]. The 
obtained upper limits scale linearly with the beaming 
factor (since we only expect to see sources for which the 
beam points towards us). 

We calculate population upper limits for a range of 
GW isotropic emission Eqw and neutrino emission n UEN . 
The results are shown in Figure 1, both for initial and 
advanced detectors. 

To interpret the science reach of the expected 
GW+HEN population upper limits described above, we 
consider the source parameters from some of the emis- 
sion models, as discussed in section II. In Figure 1 hori- 
zontal lines indicate the neutrino rate predictions of the 
Waxman-Bahcall emission model [20] , as well as the emis- 
sion model for mildly relativistic jets by Horiuchi and 
Ando [51]. The population upper limit estimates for 
these two models specifically, as functions of Eqw down 
to E GW = 10 -4 M c 2 , are shown in Figure 2. For sources 
of weaker GWs than 10 -4 M c 2 as predicted by some 
CCSN simulations [35, 41, 42], observations will focus on 
galactic sources. 
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FIG. 2. Expected GW+HEN source population upper limits 
for anticipated HEN emission from two emission models, as 
functions of isotropic-equivalent GW emission energy Eq^. 
Results are shown both for measurements with the initial 
LIGO- Virgo detectors and the IceCube-22 detector (dashed 
line), as well as for the advanced LIGO- Virgo detectors and 
the IceCube-86 detector (solid line), both with one year of 
coincident measurement time. For comparison, the galactic 
supernova rate is shown (dotted line). This Figure shows a 
subset of the results shown in Figure 1. 



IV. DATA 

In this section we describe the output of the GW and 
HEN detectors, as well as the astrophysical source dis- 
tribution. We derive the quantities that will be used in 
joint GW+HEN analyses. 



A. Gravitational- wave Data 

Interferometric GW detectors measure a passing GW 
by monitoring the distance between test masses using 
laser interferometry [73] . The detectors considered in the 
present analysis are Michelson interferometers, in which a 
laser beam is divided into two perpendicular laser beams 
directed along the two arms of the interferometer. The 
two arms, with roughly equal length L, contain resonant 
Fabry-Perot optical cavities with partially transmitting 
input mirrors and highly reflective end mirrors [1]. A 
passing GW changes the phase of the laser light between 
the entering and exiting beams with 180° phase shift be- 
tween the two perpendicular arms. The induced phase 
shift is measured through the interference of the outcom- 
ing beams from the two arms, and is used to calculate 
the so-called differential arm length AL = L\ — L 2 with 
Li and L 2 being the lengths of the two arms. The quan- 
tity h(t) = AL/L is called the GW strain, and is used to 
define the amplitude of a GW. 



The measured GW strain h(t) is denned by the passing 
GWas 

h(t) = F+h+(t) + F x h x (t) (3) 

where + ("plus") and x ("cross") are the two polariza- 
tions of the GW, at 45° angles from each other. The 
coefficients F+ and F x are the so-called antenna factors 
that depend on the direction of the incoming GW rel- 
ative to the orientation of the detector, as well as the 
polarization of the GW (see, e.g. [74]). h+{t) and h x (t) 
are the amplitudes of the GW in the two polarizations. 

GW search algorithms are designed to detect and ex- 
tract information about a GW signal from a stream of 
strain data from a set of GW detectors (e.g. [75-77]). 
One can think of a generic search algorithm as a Glf 
radiometer, outputting the excess GW energy measured 
by a network of detectors, as a function of time t and sky 
location xl. These data analysis algorithms usually out- 
put so-called GW triggers, potential GW signals whose 
significance exceed a given threshold. A GW trigger's 
significance is characterized by a suitable test statistic 
(see e.g. [75, 76]). GW triggers can have additional pa- 
rameters, such as time of arrival, amplitude or waveform. 
Data from a network of GW detectors can also be used 
to recover directional information (e.g., [78]). 

For the purposes of the joint GW+HEN analysis, we 
consider short transient events. The duration of a tran- 
sient GW event is expected to be much shorter than the 
coincidence time window [69] of GW and HEN events 
(the window in which all GW and HEN signals arrive). 

To obtain the background distribution of GW triggers, 
we time-shift data from the different GW detectors com- 
pared to each other such that no astrophysical signal can 
appear in more than one detector at a time. Background 
triggers can be generated in such a way for many different 
time shifts. 

Let us assume that we have a GW search algorithm 
that identifies a set of GW triggers, for each trigger cal- 
culating a test statistic TS and a skymap (point spread 
function) 7" GW (^). The point spread function gives the 
probability distribution of source direction, given that 
the GW event is of astrophysical origin. To calculate the 
significance of a joint event, we need to take into account 
TS as well as J GW (^). The background distribution of 
TS can be obtained from time-shifted data. The distri- 
bution of TS for the case of a signal present, however, is 
not available, as it would greatly depend on the proper- 
ties and direction of the signal. Therefore, we take into 
account TS in the joint significance by calculating its p- 
value, given the background distribution. 

Let FAR^ be the false alarm rate of GW event z, corre- 
sponding to the rate of GW events with TS>TSi (average 
number of events over unit time). For TS^ we assign a 
p- value of 

p^ = l-Pois(0;T.FAR;), 

where Pois(fc, A) is the Poisson probability of k outcome 
with A average, and T is the coincidence time window (see 
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section V A) . FAR^ is calculated from the distribution of 
background events. 

For the skymap J" GW (^) both the background and sig- 
nal distributions are available. We therefore take this 
information into account in the form of a likelihood ra- 
tio. Here our null hypothesis is that the GW event is a 
random fluctuation from the background, i.e. it has no 
preferred direction. We therefore approximate the back- 
ground likelihood Bqw to be a flat distribution over the 
whole sky, i.e. 

B« = l- (4) 

Our alternative hypothesis is that GW event i came from 
an astrophysical source at direction x^. The signal like- 
lihood <SgvI t> e the calculated skymap, i.e. 

B. High-energy Neutrino Data 

HENs traveling through the Earth interact with the 
surrounding matter with a small interaction cross section. 
In charged-current interactions, most of the neutrino's 
energy is carried away by a single high-energy electron, 
muon, or tau particle, which will emit Cherenkov radi- 
ation as it travels through the detector medium (water 
or ice). For neutrino astronomy, the high-energy muons 
are generally the most useful: they neither lose energy as 
rapidly as electrons nor decay as rapidly as taus, and 
therefore can have paths many kilometers long. The 
Cherenkov light emitted along this path can be detected 
and used to measure the direction and energy of the muon 
and infer similar information about the primary neutrino. 
HEN detectors contain large numbers of optical sensors 
placed along vertical wires (strings). These optical sen- 
sors detect the Cherenkov photons emitted by muons. 
z 

The direction of HENs can be reconstructed using the 
arrival time of Cherenkov photons at different optical sen- 
sors, with a precision of ~ 0.5° — 1° (depending on en- 
ergy) for IceCube [79], or less than 0.3° for Antares [80]. 
Direction reconstruction is also one of the major tools 
in background rejection. So-called atmospheric muons, 
created by cosmic rays interacting with particles in the 
atmosphere over the detector, are the dominant back- 
ground. To suppress these events, searches for neutri- 
nos are principally performed using up-going events, i.e. 
those that have traveled through Earth and therefore can 
be attributed only to a neutrino. The vast majority of 
these up-going neutrinos are themselves the result of cos- 
mic ray interactions on the other side of the Earth. These 
are the so-called atmospheric neutrinos, and in general 
constitute an irreducible background in searches for as- 
trophysical neutrinos from space. 

Many sources of astrophysical neutrinos are expected 
to exhibit a harder energy spectrum (typically E~ 2 ) com- 
pared with the soft spectrum of atmospheric neutrinos 



(~ E~ 3 - 7 ). In such cases, reconstructed energy can play 
a role in further separating signal from background. The 
measured energy of the muon (related to the amount of 
light detected) becomes an estimated lower bound for 
the neutrino primary energy, since the muon may have 
propagated an unknown distance before reaching the de- 
tector. Probability distribution functions for the muon 
energies expected from signal (of different source spectra) 
and background can then be used in likelihood analyses 
to enhance sensitivity for hard source spectra while re- 
taining sensitivity to softer spectra [81, 82]. 

A given reconstructed neutrino event i consists of a 
time of arrival t„\ reconstructed direction x%, direc- 
tional uncertainty o~i and reconstructed neutrino energy 
Ei. The neutrino point spread function, i.e. the prob- 
ability distribution of the neutrino source direction, is 
defined as 

M^i) = ^-e~^T^, ( 6 ) 

where x~t is the true sky location of the source. The 
HEN point spread function is incorporated in the joint 
GW+HEN significance in the form of a likelihood ratio. 
Our null hypothesis is that HEN event i is a detected 
atmospheric neutrino, having no preferred direction. We 
approximate the background likelihood B^ to be a flat 
distribution over one hemisphere (since a neutrino detec- 
tor is only sensitive to roughly half the sky), i.e. 



Our alternative hypothesis is that the neutrino came from 
an astrophysical source at direction x~l. The signal like- 
lihood S$ will be the point spread function, i.e. 

5»(5j) = Jv(xt). (8) 

For reconstructed neutrino energy the background 
distribution is known from the detected (background) 
neutrinos. The distribution of Ei for astrophysical sig- 
nals, however, depends on the source emission model, 
therefore we treat it as unknown. We therefore take into 
account Ei in comparison with its background distribu- 
tion by calculating its p- value Phen? defined as the fraction 
of background neutrinos having energy E > E^. 

C. Neutrino Clustering 

After obtaining the properties of individual neutrino 
events, we consider the possibility that multiple neutrinos 
are detected from the same astrophysical source. A set of 
neutrinos can potentially originate from the same source 
only if they are spatially and temporally coincident. 

As it is described in section V A below, assuming that 
a GW+HEN source emits neutrinos within a time in- 
terval At u (denoted by t^ — ti ^ in section VA), we 
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consider a set of neutrinos to be temporally coincident if 
all neutrinos arrive within a At v time interval. 

For spatial coincidence, we require that all neutrinos 
have a common direction of origin from where each neu- 
trino can originate with probability above a threshold 
Pmin = 0-05 (i.e. if the probability that a neutrino 
came from a direction farther from its reconstructed di- 
rection than the common direction is < P m in)- This 
probability threshold corresponds to an angular differ- 
ence threshold Alzf™ a;E between the common direction 

and the neutrino direction (AH?™^ = ais/2 ln(l/P m ^ n ). 
The average angular distance threshold for IceCube-22 
neutrinos is ~ 2.5° for P m in = 0.05). Neutrinos i and 
j will be spatially coincident if their angular distance is 
< A~af™? x + A1?™J X . 

To take into account more than one neutrino in the 
joint GW+HEN significance, one needs to account for 
the probability of detecting a certain number of neutri- 
nos from the background. The distribution of the number 
of neutrinos expected for astrophysical signals is model- 
dependent due to the non-homogenous source distribu- 
tion, and is in general unknown. However, the probabil- 
ity of background neutrinos to occur in the same time 
window can be calculated. Given that we have at least 
one detected neutrino, the probability of detecting N neu- 
trinos in the allowed time window is Pois(7V — 1; f v At v ), 
where f v is the background neutrino rate (its typical 
value is - 10~ 4 Hz for IceCube-22 and - 10~ 2 - 5 Hz for 
the completed IceCube). We calculate the p- value for 
one neutrino to be in a time-window with N or more 
neutrinos to be 

N-2 

Pcluster(N) = 1 - P ° iS (^ fvAt v ). (9) 

i=0 

This p-value is additionally taken into account to the 
other p- values in the joint test statistic. Note that, 
if one has only one detected neutrino in the cluster, 

PclusteriX) 1" 

Note that the decision of whether to treat coincident 
neutrinos as a cluster or as individual events is made by 
the analysis based on which combination yields the high- 
est significance (this decision process can proceed itera- 
tively on the remaining neutrinos until all are accounted 
for). 

D. Astrophysical Source Distribution - the Galaxy 
Catalog 

The distribution of astrophysical GW+HEN sources 
at detectable distances is not uniform. This can be 
used in a joint search algorithm to increase sensitivity. 
One can weigh event candidates based on the expected 
source density in their direction. The density of pro- 
posed GW+HEN sources can be connected to the blue 
luminosity of galaxies [83, 84], while source density can 
also depend on, e.g., the galaxy type [85, 86]. We take 
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FIG. 3. Probability that a random sky direction falls by 
chance within a given angular distance of at least one of the 
nearby galaxies listed in the GWGC. 

the blue-luminosity distribution of galaxies up to 40 Mpc 
from the Gravitational Wave Galaxy Catalogue (GWGC) 
[19]. 

Let the astrophysical GW+HEN source density be 
p(r, of), where r and ~£ = (0, 6) are the source distance 
and direction on the sky, respectively. The probability 
distribution of astrophysical neutrinos as a function of 
direction is proportional to the number of sources in the 
given direction, weighted with the distance of the sources 
to the —2 nd power (which cancels out in the volume in- 
tegral): 

Fgaii^s) = -rf- / p(r,^t)dr (10) 

Jo 

where N v is a normalization factor, xl is the source di- 
rection and Dhorizon is an expectation-motivated cut- 
off distance (see, e.g., [19, 87, 88]). For HEN searches, 
Dhorizon can be chosen to be very large, however for joint 
GW+HEN searches Dh orizon will be chosen to be the 
cutoff distance related to GW detection (see below). For 
searches using initial GW detectors, a reasonable choice 
can be Dh orizon = 40 Mpc. Given the detector sensi- 
tivities and typical source strengths, sources of interest 
father than this distance are unlikely to have measurable 
effect. In the following we will refer to this distribution 
as the weighted galaxy distribution. 

To take into account the galaxy distribution in the joint 
analysis, we consider our null hypothesis to be that the 
joint signal is a random coincidence from the background, 
i.e. it has no directional preference. This results in a 
background likelihood of 



where we take into account that a joint event can only 
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come from half of the sky due to the directional sensitiv- 
ity of neutrino observatories. The alternative hypothesis 
is that the joint event came from an astrophysical source 
at direction xl. The corresponding signal likelihood is 

Sfaltfs) = Fgal&). (12) 

The information on the distribution of galaxies is ac- 
curate for directions outside the galactic plane. Within 
the plane, the large density of galactic stars makes it 
more difficult to detect galaxies in these directions. This 
incompleteness needs to be taken into account in our per- 
ceived source distribution. Another complication is that 
nearby galaxies can be smeared to a finite area of the 
sky. We ensure that no source is missed due to these 
incompletenesses by performing a complementary search 
with no galactic weighing. Such a search, while being 
significantly less sensitive than a search that takes into 
account the galaxy catalog, is capable of detecting strong 
sources that are not aligned with the galaxy catalog. For 
this case, the galactic likelihood ratio is uniformly taken 
to be unity. 

To illustrate the capabilities of using the galaxy cat- 
alog in rejecting false GW+HEN coincidences, we cal- 
culated the probability that a random sky direction is 
within a given angular distance from at least one galaxy 
in the galaxy catalog. This is a simpler and less sensi- 
tive way of utilizing information on galaxy locations than 
used in the method (which includes the blue-luminosity 
weight), but it already shows the usefulness of this addi- 
tional information in background rejection. The results 
are shown in Figure 3. The probability of the coinci- 
dence of a random sky direction and at least one galaxy 
is evaluated for angular distances ranging from 0.1 to 10 
degrees and considering galaxies with four different hori- 
zon distances from the observer. For these curves, it is 
possible to estimate the probability that a background 
neutrino be falsely associated with a host galaxy. We see 
that for an horizon distance of 50 Mpc (which is larger 
than the 40 Mpc used in the present analysis), there is 
one chance in two to get a false association if the position 
uncertainty is of order of one degree. The efficiency of 
the galaxy catalog to discard background neutrino trig- 
gers is directly connected to this probability (since back- 
ground neutrinos coming from directions where there is 
no galaxy can be discarded). This result indicates that, 
within the 40 Mpc distance horizon, most background 
GW+HEN events will be farther from any galaxy than 
the typical angular uncertainty of HEN direction recon- 
struction, demonstrating the benefit of using the galaxy 
catalog, even for this simple model. 

V. JOINT GW+HEN ANALYSIS 

This section describes the joint analysis method for 
the search for GW+HEN signals. The joint analysis is 
described with a flow diagram in Fig. 4. For easier nav- 
igation the different steps in the flow diagram include 



references to the sections and figures where they are de- 
scribed in detail. 

While the presented search example focuses on GWs 
and HENs, we note that it is straightforward to use the 
method with other messengers as well. The method also 
naturally incorporates externally triggered searches (see 
e.g. [89, 90]), where at least one messenger confirms the 
presence of an astrophysical signal. A confirmed signal 
can define either or both the time window and source lo- 
cation (or point spread function) which restricts the pa- 
rameter space of the multimessenger search a-priory. For 
such cases the interesting scientific question is whether 
additional information is present in other messengers, 
and if it is then what can one infer about the source 
by the total available information. We also note that 
the search can be analogously used for much lower en- 
ergy neutrinos. For instance because the photomultipier 
tube (PMT) dark noise rate is particularly low in ice, 
the IceCube detector has sensitivity to sudden fluxes of 
MeV neutrinos which lead to collective rise in the PMT 
rates. Nearby supernovae up to 50 kpc are expected to 
be detected this way. While the MeV neutrino signal 
does not provide any directional information, it can be 
readily naturally incorporated in the present joint analy- 
sis by using its time of arrival and significance (i.e. flux). 
The lack of directional information can be taken into ac- 
count as a uniform sky distribution. Further, similarly 
to the blue-luminosity-weighted galaxy distribution, a- 
priori source source distribution can be used for nearby 
sources as well, for example in the form of the matter dis- 
tribution within the Milky Way. Galactic sources behind 
the center of the Milky Way can be especially interest- 
ing for multimessenger searches since they are difficult to 
observe electromagnetically. 



A. Coincidence Time Window 

The maximum time difference between the arrivals of 
the observed GW trigger and HEN events is one of the 
key parameters of the joint GW+HEN search algorithm 
[69]. A too small time window might exclude some po- 
tential sources, while a too large time window would un- 
necessarily increase the false alarm rate and the compu- 
tational cost. 

Here, we adopt a conservative arrival time difference 
of ±500 s derived for GRBs by Baret et al [69]. Given 
a neutrino event, this allows for ±500 s coincidence time 
window for a GW trigger. Multiple neutrino events and 
a GW trigger are considered temporally coincident if the 
greatest time difference between any two of these neu- 
trinos, or any neutrino and the GW trigger, is less than 
500 s. 

We consider only one GW transient (trigger) per as- 
trophysical GW±HEN source (we choose the GW trigger 
that gives the maximum joint significance; see below). 
Besides determining temporal coincidence, we apply no 
additional weight based on the arrival times of the HEN 
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FIG. 4. Flow diagram of the joint GW+HEN search algorithm. Steps include references to the sections and/or Figures in 
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events and GW trigger (while the flux of neutrinos is 
probably time dependent, the uniform weight reflects our 
lack of information about this time dependence; see, e.g, 
[91]). 



B. Joint GW+HEN Significance 

The joint significance combines the significances of the 
GW, HEN, and galaxy distribution components. For the 
directional distribution of these components, there exists 
a signal hypothesis, and therefore they are assigned a 
likelihood ratio. The other components of the GW and 
HEN events (e.g. energy) have no model-independent 



signal hypotheses, therefore they are assigned a p- value. 
These two types of information are combined into one 
joint significance measure. 

First, we combine the likelihood ratios of the direc- 
tional components to obtain a significance measure, i.e. 
p-value. The joint likelihood ratio £(xl) is defined as 



s^s)s^ s )U {j} si J \^s) 



(13) 



where {j} is the set of neutrinos within GW+HEN trig- 
ger i. Note that the joint likelihood ratio, as it combines 
the directional distributions, is defined as a function of 
direction. Example directional distributions are shown in 
Fig. 5. Since we are mainly interested in the significance 
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of the signal being of astrophysical origin, we treat the 
direction as a nuisance parameter and marginalize over 
it. Since the background likelihoods are uniform over the 
sky, the marginal likelihood ratio is 



(14) 



The background distribution of £, P B g(£), can be ob- 
tained from time scrambled data (see section VC). Using 
this distribution, the p-value p s k y of the directional part 
of the joint event can be calculated as 



(0 

Psky 



rOO 

/ P BG (C')dC. 
Jew 



(15) 



We follow Fisher's method in combining the p-values 
into one joint test statistic: 



Xf = -2 In 



PakyPolp Cluster -(N) J\p^ 
{j} 



(16) 



To ensure that possible correlations between the differ- 
ent p-values do not affect the outcome, we calculate the 
significance of X 2 by calculating its p-value from its back- 
ground distribution P BG (X 2 ) in time-scrambled data (see 
section V C): 



Pgwhen — / Pbg(X 2 )dX z 
Jx 2 



(17) 



Note that, for a given joint event, we consider the single 
HEN or cluster of HENs with only one GW trigger at a 
time. If there are more GW triggers coincident with a 
given HEN trigger, we treat them as separate joint events 
(combined with the same neutrino). 

To take into account possible sources outside of galax- 
ies and in the Milky Way (or in galaxies not included 
in the used galaxy catalog), we perform an additional 
search without using the galaxy distribution. This is 
done simply by taking T ga \ to be unity over the whole 
sky. This additional search without the galaxy catalog is 
represented by the on/off switch of the galaxy informa- 
tion in the flow diagram in Figure 4. 



C. Background Trigger Generation 

In order to calculate the significance of one or more 
joint signal candidates, we compare their test statistic X 2 
to the test statistic distribution of background coincident 
triggers. The steps of the generation of the background 
distributions are described below, and are also shown in 
Figure 6. 

For a set of GW and HEN detectors, we apply time 
shifts for background event generation. For GW de- 
tectors the time shifts between any two detector data 
streams are selected to be much greater than the maxi- 
mum possible time shift for an astrophysical signal. For 



neutrino detectors the time of arrivals are scrambled be- 
tween neutrinos, keeping each event's local coordinates 
(#, (j)) and energy fixed during the scrambling. This pro- 
cedure reproduces fairly well the distribution of back- 
ground neutrino parameters, and also preserves the geo- 
metric acceptances of the GW and HEN detectors which 
are fixed with respect to each other. 

Using the background data described above, one can 
calculate the test-statistic distribution of background 
triggers similarly to how it is done for real data with 
no time shifts (see Eq. (17)). 



D. Individual Detection 

The loudest GW+HEN event in real data will be con- 
sidered a joint detection if its probability of arising from 
the background during a one year long measurement pe- 
riod is less than 2.87 x 10 _T (one-sided 5a). 

We consider the joint GW+HEN to have discovery po- 
tential for a given GW+HEN flux measured at Earth 
from an astrophysical source if such a signal has 50% 
probability of resulting in a joint detection (as defined 
above) . 

We define the median upper limit of the joint 
GW+HEN search to be the joint GW+HEN flux at 
Earth, in terms of h rss and (n), for which X 2 is greater, 
with 90% confidence, than the median of the loudest 
events among sets of N randomly generated joint back- 
ground events. 

To evaluate the statistical significance of a single 
GW+HEN event with X 2 , we compare it to the back- 
ground X 2 distribution P BG (X 2 ) as described in section 
VC. The statistical significance of a X 2 value is given by 
its p-value (Eq. 17). 



E. Statistical Detection of Multiple Sources 

Besides detecting a single astrophysical GW+HEN 
joint event, one can try to indicate the presence of multi- 
ple astrophysical joint events statistically. For such sta- 
tistical detection, we compare the distribution P S (X 2 ) of 
the real data to the distribution P BG (X 2 ) of joint back- 
ground events. We use the realistic assumption that only 
a small fraction of the signal candidates can be due to 
actual astrophysical signals. The steps of statistical de- 
tection described below are also shown in Fig. 7. Two 
alternative statistical tests can be found, e.g. in [92]. 

If only a small fraction of the joint events in real data 
are from astrophysical sources, one has the best chance 
of detecting the presence of real astrophysical signals by 
looking at the highest X 2 values. We therefore select a 
X 2 threshold above which the real-data and background 
distributions are compared. We denote this threshold 
by X 2 . This threshold is chosen based on the back- 
ground neutrino event rate, and the estimated astrophys- 
ical neutrino event rate within the distance in which the 
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FIG. 5. Example likelihood distributions as parts of the joint likelihood ratio: weighted galaxy distribution (upper left), HEN 
directional probability distribution function (PDF) (upper right), GW PDF (lower left) and joint PDF (lower right). The joint 
PDF is the product of the other three PDFs. The scales are in arbitrary units, and the color-scale for the galaxy distribution 
is inverted. On the joint PDF plot, every galaxy for which the joint PDF is non-zero is circled for visibility. The reconstructed 
source direction with the maximum significance is circled with bold line. 
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GW+HEN search is sensitive. 

We compare the distributions of real and time-shifted 
data above X? in the following manner. Let p t be the 
p- value corresponding to Xf. We introduce the product 
p for a set of p- values which are above threshold p t . The 
value p can be written as 



] [ Pgwhen 



(18) 



Pgwhen ^>Pt 



where Pgwhen is the p-value of measurement i (see Eq. 
17). Similarly to the use of p- values for the single detec- 
tion case, we calculate the probability p p that the mea- 
sured p from real data can arise from the background: 



f 

Jo 



(19) 



where P B c(p) is the probability distribution of p on time- 
shifted data (of identical duration as real data). The 
value p p is therefore the probability that the product of 
the p- values smaller than p t from real data arose from the 
background. We use p p to characterize the significance 
of statistical detection. We claim statistical detection 
if the probability that the real-data p p arose from the 
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background during one year of measurement is less than 
2.87 x 1(T 7 (one sided 5a). 



F. Simulation of Astrophysical Signals 

We use simulated astrophysical signals to character- 
ize the sensitivity of the GW+HEN search algorithm. 
The simulation is designed such that the results are scal- 
able for different GW and HEN emission fluxes, there- 
fore our results can be used to constrain the parameter 
space of GW and HEN emissions. We simulate standard 
GW+HEN sources with identical intrinsic GW+HEN 
emission (energy and spectrum). Upper limits for such 
standardized source populations are conservative esti- 
mates, taking the average emission, compared to taking 
a distribution of emission energies. 

For the simulations we assume that zero or one HEN 
event is detected from each source. This is the likely 
situation for the part of the HEN parameter space that 
is not constrainable by a km 3 HEN detector alone. 

The simulation of an astrophysical joint event consists 
of the following steps. 

1. For a given direction, generate a simulated astrophys- 
ical HEN event coming from the source direction. We 
use Monte Carlo simulations to generate a random 
reconstructed energy and directional uncertainty for 



such a neutrino, using a source neutrino energy spec- 
trum. We then generate a reconstructed source direc- 
tion for the neutrino, drawn from a 2D Gaussian distri- 
bution centered around the real source direction, with 
standard deviation obtained from the Monte Carlo 
simulation. 

2. Generate simulated astrophysical GW event coming 
from the source direction. We inject a GW signal with 
a given amplitude and waveform into real GW data 
into each GW detector used in the analysis, taking into 
account the direction of the source. The amplitude of 
the injected GW is chosen from a range that covers 
the amplitude region of interest for the joint search. 



G. Population Estimation 

We define population upper limit for the joint emitters 
of GWs and HENs as the lowest population which would 
produce - with > 90% probability, a joint event with 
higher significance than the loudest GW+HEN event in 
real data. 

We obtain population estimates by calculating the 
probability of detection rate for every galaxy. We adopt 
the following notation: for an astrophysical source in 
galaxy z, the galaxy has blue luminosity L^) and is at 
a distance T{. Further, let be the blue luminosity 

of the Milky Way, R the source rate [per year per MWE 
galaxy], T m the measurement duration, and ft the neu- 
trino beaming factor of the source. The calculation is the 
following. 

Given an astrophysical source in galaxy z, the proba- 
bility that at least one HEN will be detected from this 
source is [93] 

p(n > l|r,n H E N ) = l-F poiss (0|(10 Mpc/r) 2 

^HEN 

. (20) 

where F poiss is the Poisson cumulative distribution func- 
tion, and n is the number of detected neutrinos from the 
source. Therefore for galaxy i the average number Ni of 
sources with at least one detected neutrinos during the 
measurement will be 

N i (R,T m )=p(n > l^un^-R/h-T-L^/Lf™. (21) 

The population upper limit is obtained from Ni(R, T m ) 
by requiring the total number of detected neutrinos 
within £) GWHEN to be 2.3 during a one year long mea- 
surement. This is done by summing Ni(R,T m ) over all 
galaxies on the hemisphere in which the neutrino detec- 
tor is sensitive. For IceCube this is 5{ > where 5{ is the 
declination of galaxy i. The population upper limit will 
be 
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where e GW (ri) is the detection efficiency of the GW de- 
tector network at distance r$ [72]. 



H. Estimated Sensitivity 

In Section III we estimated the science reach of the 
baseline multimessenger analysis following the calcula- 
tions of Bartos et al. [93]. In estimating the science 
reach the single required parameter from the search al- 
gorithm was its exclusion distance D™% EN . We approxi- 
mated this parameter for the analysis by assuming that 
it will be comparable to the horizon distance of exter- 
nally triggered GW searches. In the externally trig- 
gered GW search of Abbott et al. [16], the authors give 
the median horizon distance for single GW events co- 
incident with electromagnetically observed GRBs. This 
horizon distance is related to the expected loudest GW 
background event from a given direction within a given 
time window of ~ 100 s. The approximation of D^™ N 
for the baseline multimessenger search is reasonable if, 
given the number of HENs, the GW triggers' significance 
threshold is chosen such that the expected number of 
spatially and temporally coincident GW+HEN events is 
< 1. For one coincident event the remaining difference 
between the GW-GRB externally triggered search and 
the GW+HEN multimessenger search is mainly due to 
the greater directional uncertainty of neutrinos (~ 1°) 
compared to the much better directional accuracy of elec- 
tromagnetic GRB measurements. This difference, how- 
ever, will not be significant, as the directional accuracy 
of GW measurements [0(10°)] is much worse that of 
the HEN directional accuracy (see e.g. [78, 94]). Fur- 
ther, the GW+HEN multimessenger analysis addition- 
ally takes into account the significance of HENs based on 
their reconstructed energy (see Section IV B) as well as 
the expected source distribution (see Section IV D). Both 
of these pieces of additional information further increase 
the sensitivity of the search, making the comparison to 
results from externally triggered GW searches conserva- 
tive. 

To estimate the validity of the approximation that one 
will have < 1 joint event in a measurement without signif- 
icant constraints on the rate of GW triggers, we take the 
example of the initial LIGO- Virgo detectors during their 
S5/VSR1 science run and the partially completed Ice- 
Cube detector in its 22-string configuration. The LIGO- 
Virgo-IceCube network ran in coincidence from May 31, 
2007 until Sept. 30, 2007 (123 days). IceCube-22, dur- 
ing its full run of 275.7 days, collected a final sample 
of 5114 neutrino candidate events [95], of which ~ 1000 
occur during the coincident livetime of the LIGO-Virgo- 
IceCube network. Considering a characteristic GW point 



spread function that spreads over 100 deg 2 , the proba- 
bility that temporally coincident GW and HEN triggers 
are also spatially coincident is ~ 0.5%. To reach on av- 
erage ~ 1 spatially and temporally coincident joint event 
during the measurement, we need on average 1 GW trig- 
ger for every 5 HEN triggers, corresponding to a GW 
trigger frequency of 17 day -1 . This limit is far from lim- 
iting the sensitivity of the search. For comparison, ~ 1 
GW trigger/day was used for electromagnetic follow-up 
observations for the latest LIGO- Virgo (S6/VSR2-3) sci- 
ence runs [96]. 

Comparing the expected population upper limit of 
the GW+HEN multimessenger search to an all-sky GW 
search, the multimessenger search is expected to give 
stricter constraints on source population if its increased 
horizon distance compensates for HEN beaming (GWs 
are only weakly beamed). An approximate comparison 
is given by the ratio of the number of sources above the 
expected loudest-event threshold in each of the searches: 

(new \3/f. ' 

where D™ % is the horizon distance of an all-sky GW 
search and fb, gw is the GW beaming factor. Taking 
^50% EN ~ 12 Mpc from externally triggered GW searches 
[16],°L>™ % « 7.8 Mpc from GW all-sky searches [93], 
fb « 14 (an observational estimate for low-luminosity 
GRBs; [67]) and fb,cw ~ 1.5 (estimated value for inspi- 
rals or accretion- type GW emission; e.g. [97]), the ratio 
of detectable GW+HEN and GW events is w 0.4, indicat- 
ing that the number of sources excluded with the joint 
search is comparable to the number of those excluded 
with GW all-sky searches. Further, the sources probed 
by a joint search are mostly complementary to sources 
probed by an all-sky GW search. As the joint analysis is 
looking farther, it can potentially see sources missed by 
GW all- sky searches. 

VI. SUMMARY 

We presented a baseline method for the joint search 
for GWs and significances of GW and HEN signals, in- 
cluding the point spread functions of event sky locations, 
as well as the blue-luminosity-weighted distribution of 
known galaxies. 

We estimated the expected science reach of a joint 
GW+HEN search based on some of the existing GW and 
HEN emission models. The expected results indicate that 
the GW+HEN search with initial and particularly with 
advanced detectors will constrain the parameter space of 
some of the existing models. 
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To interpret the science reach of the expected 
GW+HEN population upper limits, we considered HEN 
rate expectations based on the Waxman-Bahcall emission 
model [20] , as well as the emission model with mildly rel- 
ativistic jets by Ando and Beacom [51]. 

The baseline GW+HEN analysis is expected to result 
in better upper limits than independent searches if the 
increase in exclusion distance compared to an indepen- 
dent GW search is greater than ~ f^ 3 (where fa is 
the HEN beaming factor; see also [93]). Further, the 
main advantage of a multimessenger search over individ- 
ual searches is that of detection efficiency, as well as ad- 
ditional scientific information available from the source. 
Due to the non-Gaussianity of the GW data stream [1], 
high-significance events have relatively high false alarm 
rate and therefore a coincident messenger can greatly in- 
crease our confidence in a detection. Additionally, joint 
detection would help better understand the underlying 
physics of the source. For instance a HEN detected prior 
to a GW from a GRB would indicate precursor activity 
in the GRB. The time delay of the earliest HEN after 
the GW signature of a collapsar may be indicative of jet 
propagation within the stellar envelope. 
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